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Abstract. Numpy and SciPy are program libraries for the Python 
s ! scripting language, which apply to a large spectrum of numerical 

and scientific computing tasks. The Sage project provides a multi- 
^ platform software environment which enables one to use, in a uni- 

fied way, a large number of software components, including Numpy 
^--^ and Scipy, and which has Python as its command language. We 

review several examples, typical for scientific computation courses, 

and their solution using these tools in the Sage environment. 

<^ Keywords. Numerical methods teaching, Python language, Sage, 

computer algebra systems 

2010 Mathematics Subject Classification. 65-01 

-I— > 

a 

1. Introduction 

Python is a popular multipurpose programming language. When 
augmented with special libraries/modules, it is suitable also for sci- 
entitle computing and visualization. These modules make scientific 
CN computing with Python similar to what commercial packages such as 

^ MATLAB, Mathematica and Maple offer. Special Python libraries 

like NumPy, SciPy and CVXOPT allow fast machine precision floating 
point computations and provide subroutines to all basic numerical anal- 
ysis functions. As such the Python language can be used to combine 
several software components. For data visualization there are numer- 
ous Python based libraries such as Matplotlib. There is a very wide 
^ collection of external libraries. These packages/libraries are available 

either for free or under permissive open source licenses which makes 
these very attractive for university level teaching purposes. 

The mission of the Sage project is to bring a large number of program 
libraries under the principle of the GNU General Public License under 
the same umbrella. Sage offers a common interface for all its compo- 
nents and its command language is Python. Originally developed for 
the purposes of number theory and algebraic geometry, it presently 
supports more than 100 special libraries. Augmented with these, Sage 



is able to carry out both symbolic and numeric computing. Sage 



File: tmj_sage24.tex, printed: 2012-8-21, 0.29 

1 



2 



In this paper, we present examples or case studies of the usage of 
Sage to solve some problems common for typical numerical analysis 
courses. The examples are drawn from the courses of the second author 
at the University of Turku, covering the standard topics of numerical 
computing and based, to a large extent, on the standard textbooks 

[HB EdHl El 1MB M E El [tln] . 

2. General observations 

2.1. History of Sage. Sage was initially created by William Stein in 
2004-2005, using open source programs released under the GPL or a 
GPL-compatible license. The main goal of the project was to create a 
viable open source alternative to proprietary mathematical software to 
be used for research and teaching. The first release of the program was 
in February 2005. By the end of the year, Sage included Pari, GAP 
and Singular libraries as standard and worked as a common interface 
to other mathematical programs like Mathematica and Magma. [His, 

Since the beginning Sage has expanded rapidly. As of July 2012, at 
least 246 people have actively contributed code for Sage. The range of 
functionality of Sage is vast, covering mathematical topics of all kinds 
ranging from number theory and algebra to geometry and numerical 
computation. Sage is most commonly used in university research and 
teaching. On Sage's homepage there are listed more than one hundred 
academic articles, books and theses in which the program has been 
involved. 

2.2. Access to Sage. Sage can be utilized in multiple ways. The main 
ways to use Sage are: 

(1) Over the network. In this case no installation is required. 

(2) The software is downloaded and installed on the personal work 
station. 

These ways to use Sage will now be described more closely. 

(1) One of the strengths of Sage is that it can be used over the 
network without requiring any installation of the application. 
The Sage Notebook is a web browser-based graphical user in- 
terface for Sage. It allows writing and running code, displaying 
embedded two and three dimensional plots, and organizing and 
sharing data with other users. The Sage Notebook works with 
most web browsers without the need for additional add-ons or 
extensions. However, Java is required to run Jmol, the Java 
applet used in Sage to plot 3D objects. 
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(2) Sage can be also installed natively for Linux, OS X and Solaris. 
Both binaries and source code are available for download on 
Sage's homepage. In order to run Sage on the Windows operat- 
ing system the use of virtualization technology (e.g. VirtualBox 
or VMware) is required. There are three basic interfaces to ac- 
cess Sage locally: the Sage Notebook on a web browser, a text- 
based command-line interface using IPython, or as a library in 
a Python program. |PG] 

In addition to these main options, Sage can be applied in alternative 
ways. For instance, a single cell containing Sage code can be embedded 
in any webpage and evaluated using a public single cell server. There 
is also support to embed Sage code and graphics in LaTeX documents 



using Sagetex. SageTeX 



2.3. Key issues. Numerical computation has been one of the most 
central applications of the computer since its invention. In modern 
society, the significance of the speed and effectiveness of the computa- 
tional algorithms has only increased with applications in data analysis, 
information retrieval, optimization and simulation. 

Most numerical algorithms are implemented in a variety of program- 
ming languages. There are various commercial collections of software 
libraries with numerical analysis functionality, mostly written in For- 
tran, C and C++. These include, among others, the IMSL and NAG 
libraries. Computer algebra systems and applications developed spe- 
cially for numerical computing include MATLAB, S-PLUS, Mathemat- 
ica, IDL and Lab VIEW. There are also free and open source alterna- 
tives such as the GNU Scientific Library (GSL), R, Scilab, Freemat 
and Sage. 

Many computer algebra systems (e.g. Maple, Mathematica and 
MATLAB) contain an implementation of a new programming language 
specific to that system. In contrast, Sage uses Python, which is a pop- 
ular and widespread high-level programming language. The Python 
language is considered simple and easy to learn while making the use 
of more advanced programming techniques, such as object-oriented pro- 
gramming and defining new methods and data types, possible in the 
study of the mathematics. Python functions as a common user inter- 
face to Sage's nearly one hundred software packages. 

Some of the advantages of Sage in scientific programming are the 
free availability of the source code and openness of development. Most 
commercial software do not have their algorithms publicly available, 
which makes it impossible to revise and audit the functionality of the 
code. Therefore the use of the built-in functions of these programs 
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may be inadequate in some mathematical studies that are based on the 
results given by these algorithms. This is not the case in open source 
software, where the user can verify the individual implementations of 
the algorithms. 

In many situations, the Python interpreter is fast enough for com- 
mon calculations. However, sometimes considerable speed is needed in 
numerical computations. Sage supports Cython, which is a compiled 
language based on Python that supports calling C functions and declar- 
ing C types on variables and class attributes. It is used for writing fast 
Python extension modules and interfacing Python with C libraries. 
Significant parts of the NumPy and SciPy libraries included with Sage 
are written in Cython. In addition, techniques like distributed com- 
puting and parallel processing using multi-core processors or multiple 
processors are supported in Sage. [BBSE| INum] 

Sage includes the Matplotlib library for plotting two-dimensional 
objects. For three-dimensional plotting, Sage has Jmol, an open-source 
Java viewer. Additionally, the Tachyon 3D raytracer may be used for 
more detailed plotting. Sage's interact function can also bring added 
value to the study of numerical methods by introducing controllers 
that allow dynamical adjusting of the parameters in a Sage program. 
pirn iJnidl [Tl [Tic] 

2.4. Numerical tools in Sage. Sage contains several program li- 
braries suitable for numerical computing. The most substantial of 
these are NumPy, SciPy and CVXOPT, all of which are extension 
modules to the Python programming language, designed for specific 
mathematical operations. In order to use these packages in Sage, they 
must be first imported to the Sage session using the import statement. 
[ADVl EE |SclPy1 



NumPy provides support for fast mult i- dimensional arrays and nu- 
merous matrix operations. The syntax of NumPy resembles the syntax 
of MATLAB in many ways. NumPy includes subpackages for linear 
algebra, discrete Fourier transforms and random number generators, 
among others. The SciPy library is a library of scientific tools for 
Python. It uses the array object of the NumPy library as its basic 
data structure. SciPy contains various high level scientific modules 
for linear algebra, numerical integration, optimizing, image processing, 
ODE solvers and signal processing, et cetera. 

CVXOPT is a library specialized in optimization. It extends the 
built-in Python objects with dense and sparse matrix object types. 
CVXOPT contains methods for both linear and nonlinear convex op- 
timization. For statistical computing and graphics, Sage supports the 
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R environment, which can be used via the Sage Notebook. Some sta- 
tistical features are also provided in SciPy, but not as comprehensively 
as in R. jR] 

For more information on numerical computing with Sage, see |Numj . 

3. Case studies 

In this section our goal is to present examples of the use of Sage for 
numerical computing. This goal will be best achieved by giving code 
snippets or programs that carry out some typical numerical computa- 
tion task. We cover some of the main aspects of a standard first course 
in numerical computing, such as the books of [BF] ICdB| [Hi \MF\ IMol[ 

lEJEETN]. 

3.1. Newton's method. Computing the solution of some given equa- 
tion is one of the fundamental problems of numerical analysis. If a 
real-valued function is differentiable and the derivative is known, then 
Newton's method may be used to find an approximation to its roots. 

In Newton's method, an initial guess x$ for a root of the differentiable 
function / : R — > R is made. The consecutive iterative steps are defined 
by 

x k+ i =x k - { ^ , k = 0,1,2,... 
f{x k ) 

An implementation of the Newton's method algorithm is presented 
in the next code. As an example, we use the algorithm to find the 
root of the equation x 2 — 3 = 0. The function newtommethod is used 
to generate a list of the iteration points. Sage contains a preparser 
that makes it possible to use certain mathematical constructs such as 
f(x) = f that would be syntactically invalid in standard Python. In 
the program, the last iteration value is printed and the iteration steps 
are tabulated. The accuracy goal for the root 2h is reached when 
f(x n — h)f(x n + h) < . In order to avoid an infinite loop, the maxi- 
mum number of iterations is limited by the parameter maxn. 



def newton_method(f , c, maxn, h) : 
f(x) = f 
iterates = [c] 

j = 1 

while True : 

c = c - f (c) /derivative(f (x) ) (x=c) 
iterates . append(c) 
if f (c-h) *f (c+h) < or j == maxn: 
break 
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n 


x n 


fM 


f[x„-h)f[x„+h) 





2.000000 


1.0000 


1.000 


1 


1.750000 


0.062500 


0.003906 


2 


1.732143 


0-00031988 


1.005 x 10" 7 


3 


1.732051 


8.4727 x 10" 9 


-1.200 x 10" 9 



Figure 1. Output of the Newton iteration. 

j + = l 
return c, iterates 

f(x) = x~2-3 
h = 10~-5 
initial =2.0 
maxn = 10 

z, iterates = newton_method(f , initial, maxn, h/2.0) 
print "Root =", z 

html. table (( [i, c.n(digits=7) , f (c) .n(digits=5) , 

(f (c-h)*f (c+h)) .n(digits=4)] for i, c in enumerate (iterates) ) , 
header=["$n$", "$x_n$", "$f(x_n)$", "$f (x_n-h)f (x_n+h)$"] ) 

3.2. Computational Methods of Linear Algebra. Sage offers a 
good platform for practicing both symbolic and numerical linear alge- 
bra. The software packages specialized in computational linear algebra 
that are contained in Sage include LAPACK, LinBox, IML, ATLAS, 
BLAS and GSL. Both the NumPy and SciPy libraries have subpackages 
for algorithms that utilize NumPy's fast arrays for matrix operations. 
Also, the GSL library can be used via Cython. In most of the examples 
of this chapter, we use native Sage matrices that suit most needs. 

Let A be a matrix over the real double field (RDF) as defined in the 
next code. The matrix function accepts the base ring for the entries 
and the dimensions of the matrix as its parameters. We can compute 
various values associated with the matrix, such as the determinant, the 
rank and the Frobenius norm: 

Sage: A = matrix(RDF, 3, [1,3,-3, -3,7,-3, -6,6,-2]) 
Sage: A. determinant () 

-32.0 
Sage: A.rankQ 

3 

Sage: A.normO 
12.7279220614 
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The function A.inverse() returns the inverse of A, should it exist. 
Otherwise Sage informs that the matrix must be nonsingular in order 
to compute the inverse. 

Sage: A. inverse () 

-0.125 0.375 -0.375 \ 
-0.375 0.625 -0.375 
-0.75 0.75 -0.5 / 

Let b = (1,3,6). We solve the matrix equation Ax = b using the 
function solvejright. The notation A\b, specific to Sage, may also be 
used. 

Sage: b = vector ( [1 ,3,6] ) 
Sage: A.solve_right(b) 
(-1.25, -0.75, -1.5) 

In numerical linear algebra, different decompositions are employed to 
implement efficient matrix algorithms. Sage provides several decompo- 
sition methods related to solving systems of linear equations (e.g. LU, 
QR, Cholesky and singular value decomposition) and decompositions 
based on eigenvalues and related concepts (e.g. Schur decomposition, 
Jordan form). The availability of these functions depends on the base 
ring of matrix; for numerical results the use of real double field (RDF) 
or complex double field (CDF) is required. 

Let us determine the LU decomposition of the matrix A. The result 
of the function A.LU() is a triple of matrices P, L and U, so that 
PA = LU, where P is a permutation matrix, L is a lower-triangular 
matrix and U is an upper-triangular matrix. 

Sage: A.LUO 

/ / 0.0 0.0 1.0 \ / 1.0 0.0 0.0 \ / -6.0 6.0 -2.0 \ \ 

0.0 1.0 0.0 , 0.5 1.0 0.0 , 0.0 4.0 -2.0 

VV 1.0 0.0 0.0 / \ -0.1666 1.0 1.0 / \ 0.0 0.0 -1.3333 // 

According to linear algebra, the solution of the equation Ax = b for 
a n x n matrix is unique if the determinant det(A) ^ . However, the 
solution of the equation may be numerically unstable also if det(A) ^ . 
The standard way to characterize the "numerical nature" of a square 
matrix A is to use its condition number cond(A) defined as oul^m 
where om (°"m) is the largest (least) singular value of A . The singular 
value decomposition (SVD) of A yields the singular values as follows: 



A = USV T 
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where U and V are orthogonal n x n matrices and the S is a diagonal 
n x n matrix with positive entries organized in decreasing order, the 
singular values of A . 

In the next example we study the influence of the condition number 
on the accuracy of the numerical solution of a random matrix with a 
prescribed condition number. For this purpose we use a simple method 
to generate random matrices with a prescribed condition number c > 1: 
take a random square matrix A, form its SVD A = U SV T and modify 
its singular values S so that for the modified matrix S c the quotient of 
the largest and least singular value is c and then A c — U S C V T is our 
desired random matrix with cond(A c ) = c . For several values of c we 
then observe the error in the numerical solution of A c x = b and graph 
the error as a function of cond(A c ) in the loglog scale. We see that the 
condition number appears to depend on cond(A c ) almost in a linear 
way. 

from numpy import * 

from matplotlib.pyplot import * 

data = [] 

n = 20 

A = random. rand (n,n) 
U, s, V = linalg. svd(A) 
ss = zeros((n,n)) 

for p in arange(l , 16,2 . ) : 
c = 10. "p 
for j in range (n): 

ss[j, j] = s[0] - j * (s[0] - s[0] / c) / (n - 1) 
aa = dot (dot (U, ss) , V.T) 
b = dot(aa, ones(n)) 
numsol = linalg. solve (aa, b) 
d = linalg. norm(numsol-ones (n) ) 
data. append ( [c , d] ) 

data = array (data) 

x,y = data[: ,0] ,data[: ,1] 

clf() 

loglog (x, y, color =, k' , linewidth=2) 

loglog(x, y, 'o', color='k', linewidth=10) 

xlabelC Condition number of the matrix ',\ 

f ontweight= 'bold' , fontsize=14) 

ylabel (' Error ' , f ontweight='bold' , fontsize=14) 

title('Error as a function of the condition matrix', \ 
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f ontweight= 'bold' , fontsize=14) 

grid (True) 

savef ig ( " f ig . png" ) 



Error as a function of the condition matrix 




10 1 10 2 10 3 10* 10 5 10 6 10 7 10 8 10 9 10 10 10 11 10 12 10 13 10 14 10 ; 

Condition number of the matrix 



Figure 2. The output of the program used to study the 
influence of the condition number on the accuracy of the 
numerical solution of a random matrix with a prescribed 
condition number. 

3.3. Numerical integration. Numerical integration methods can prove 
useful if the integrand is known only at certain points or the an- 
tiderivate is very difficult or even impossible to find. In education, 
calculating numerical methods by hand may be useful in some cases, 
but computer programs are usually better suited in finding patterns 
and comparing different methods. In the next example, three numer- 
ical integration methods are implemented in Sage: the midpoint rule, 
the trapezoidal rule and Simpson's rule. The differences between the 
exact value of integration and the approximation are tabulated by the 
number of subintervals n (Fig|3|. 
f(x) = x~2 
a = 0.0 
b = 2.0 
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table = [] 

exact = integrate (f (x) , x, a, b) 

for n in [4, 10, 20, 50, 100]: 
h = (b-a)/n 

midpoint = sum( [f (a+(i+l/2) *h) *h for i in range (n)]) 
trapezoid = h/2*(f(a) + 2*sum( [f (a+i*h) for i in ranged, n)]) 
+ f (b)) 

simpson = h/3*(f(a) + sum( [4*f (a+i*h) for i in ranged, n, 2)] ) 

+ sum( [2*f (a+i*h) for i in range (2,n,2)]) + f(b)) 
table . append ( [n, h.n(digits=2) , (midpoint-exact) .n(digits=6) , 
(trapezoid-exact) .n(digits=6) , (simpson-exact) .n(digits=6)] ) 

html. table (table, header=["n", "h" , "Midpoint rule", 
"Trapezoidal rule", "Simpson's rule"]) 



n 


h 




Trapezoidal rule 


1 S3 t mn G;r*ir\ ' s tt i "1 a 1 






4 


0.50 


-0.0416666666666665 


0.0833333333333335 


0.000000000000000 


10 


0.20 


-0.0066666666666663% 


0.0133333333333341 


0.000000000000000 


20 


0.10 


-0.00166666666666604 


0.00333333333333385 


0.000000000000000 


50 


0.040 


-0.000266666666666193 


0.000533333333333275 


8.88178419700125 x lO" 16 


100 


0.020 


-0.0000666666666662152 


0.000133333333333763 


4.44089209850063 x 10" 16 



Figure 3. The table shows the difference between the 
exact value of the integral and the approximation using 
various rules. 



There are also built-in methods for numerical integration in Sage. 
For instance, it is possible to automatically produce piecewise-defined 
line functions defined by the trapezoidal rule or the midpoint rule. 
These functions can be used to visualize different geometric interpre- 
tations of the numerical integration methods. In the next example, 
midpoint rule is used to calculate an approximation for the definite 
integral of the function f(x) = x 2 — 5x + 10 over the interval [0, 10] 
using six subintervals (Fig [4]) . 

f(x) = x~2-5*x+10 

f = Piecewise([[(0,10) , f]]) 

g = f .riemann_sum(6, mode="midpoint") 

F = f .plot(color="blue") 

R = add([line([[a,0] , [a,f (x=a)] , [b,f (x=b)] , [b,0]] , color="red") 

for (a,b), f in g.list()]) 
show(F+R) 
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2 4 6 8 10 



Figure 4. The geometric interpretation of the midpoint 
rule is visualized using Sage's built-in functions for nu- 
merical analysis and plotting. 



Sage's numerical integration method numericalJntegral utilizes the 
adaptive Gauss-Kronrod method available in the GSL (GNU Scientific 
Library) library. More methods for numerical integration are available 
in SciPy's sub-packages. 

3.4. Multidimensional Newton's method. Newton's iteration for 
solving a system of equations f(x) = in R n consists of fixing a suitable 
initial value x$ and recursively defining 

x k +i = x k - J f (x k )~ l f(x k ) ,k = 0,1,2, 

Consider next the case n = 3 and the function 

(3x — cos (xix 2 ) — \ 
x 2 - 81(xi + 0.1) 2 + sinx 2 + 10.6 
e~ x ° Xl + 20x 2 + 

where x = (xq,xi,X2)- In this program the Jacobian matrix J/(x) is 
computed symbolically and its inverse numerically. As a result, the 
program produces a table of the iteration steps and an interactive 3D 
plot that shows the steps in a coordinate system. 

x0, xl, x2 = var('x0 xl x2') 

fl(x0, xl, x2) = 3*x0 - cos(xl*x2) - (1/2) 

f2(x0, xl, x2) = xCT2 - 81*(xl + 0.1)~2 + sin(x2) + 10.6 
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f3(x0, xl, x2) = e"(-xO*xl) + 20*x2 + (10*pi - 3)/3 

f(xO, xl, x2) = (f I(x0,xl,x2) , f2(x0,xl,x2) , f 3(x0,xl ,x2) ) 

j = jacobian(f, [x0,xl,x2]) 

x = vector ([3.0, 4.0, 5.0]) # Initial values 

data = [[0, x, n(norm(f (x [0] , x[l], x[2])), digits=4)]] 

for i in ranged, 8): 

x = vector ((n(d) for d in x - j(x0=x[0], xl=x[l] , 
x2=x[2]) .inverse ()*f(x[0] , x[l] , x[2]))) 
data. append ( [i , x, norm(f (x [0] , x[l], x[2]))]) 

# HTML Table 

html. table ( [(data [i] [0] , data[i] [1] .n(digits=10) , 

n(data[i] [2] , digits=4)) for i in range(0,8)] , 
header = ["$i$", "$(x_0,x_l ,x_2)$" , "$norm(f ) $"] ) 

# 3D Picture 

1 = line3d([d[l] for d in data], thickness=5) 
p = point3d(data[-l] [1] , size=15, color="red") 
show(l + p) 



D 


(.v , x^x?) 


normif) 





(3.000000000, 4.000000000, 5..000000000) 


1347. 


i 


(9.950035168, 2.038510442, - 0.4735923501) 


262.9 


2 


(0.5245370576, 0.7428920309,- 0.4735987782) 


47.13 


3 


(0.4949073626, 0.3972557032,-0.5143475312) 


9.675 


4 


(0.4971709310, 0.2771459340,-0.5170792775) 


1.169 


5 


(0.4970474904, 0.2580112360,- 0.5175788308) 


0.02966 


6 


(0.4970438097, 0.2574996297,- 0.5175919449) 


0.0001239 


7 


(0.4970438070, 0.2574992635,- 0.5175919542) 


0.0001221 



Figure 5. The program produces a table showing the 
iteration steps of the Newton's method. 



3.5. Nonlinear fitting of multiparameter functions. Given the 
data (xj,Vj),j = 1, . . . ,m, we wish to fit y — f(x, A) into a "model", 
where A = (Ai, A p ) by minimizing the object function 
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2 1 




Figure 6. An interactive 3D plot shows the iteration 
steps in a coordinate system. The plot is made with the 
Jmol application integrated in Sage. 



*(A) = Z>i-/(*i,A)) 2 - 

3=1 

The minimization may encounter the usual difficulties: the minimum 
need not be unique and there may be several local minima, each of 
which could cause the algorithm to stop prematurely. In the next 
example the function minimize uses the Nelder-Mead Method from 
the scipy. optimize package. 

Let A = (Ai, A2, A3). Consider the model function 

A) = Axe"* + \ 2 e-^ x . 

The data points used in this example are generated randomly by 
deviating the values of the model function. 
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from numpy import * 

def f model (lam, x) : 

return lam[0] *exp(-x) + lam [1] *exp(-lam [2] *x) 

def fobj (lam, xdata, ydata) : 

return linalg. norm (f model (lam, xdata) - ydata) 

xdata = arange(0, 1.15, 0.05) 
lam = [0.2, 1.5, 0.7] 
y = f model (lam, xdata) 

# The generation of the data points 
ydata = y*(0. 97+0. 05*random. rand (y. size)) 

# Initial values 
lamO = [1, 1, 1] 

yO = fobjQamO, xdata, ydata) 

# The minimization of the object function 

lam = minimize (fobj , lamO, args=(xdata, ydata), algorithm=' simplex' ) 
yfinal = fobj (lam, xdata, ydata) 

# Plot of the datapoints and the model function 

fit = plot (f model (lam, x) , (x, 0, 1.5), legend_label="Fitted curve") 
datapoints = list_plot (zip(xdata, ydata), size=20, 

legend_label="Data points") 

html ("\n\n$\\text{ Object function values: start = %s, final = °/ s}$\n" 

%( n(y0, digits=5) , n(yfinal, digits=5))) 
show(fit + datapoints, figsize=5, gridlines=True , 

axes_labels=( "xdata" , "ydata") ) 



3.6. Polynomial Approximation. The problem of finding (n — 1 )th 
order polynomial approximation for a function g on the interval [r 1? r 2 ] 
leads to the minimization of the expression 

(g(x) - c kX n ~ k Y dx 
i k=i 

with respect to the parameter vector (ci, ...,c„) . In order to find the 
optimal value of the parameter vector, we consider the critical points 
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Optimization terminated successfully. 

Current function value: 0.0682 68 

Iterations : 23S 

Function evaluations: 42 




Figure 7. The algorithm used in the program returns 
a report on the success of the optimization. The plot 
shows the data points and the model function in the same 
coordinate system. 



where gradient vanishes i.e. the points where 

df 

dci 

For the purpose of illustration, consider the case when g(x) 
n = 2, 3, 4. The equations |4 — l ea d to the requirement 



e x and 



E 

fc=i 



r 



2n-k-j+l 



2n - k - j + 1 



r 2 /-j-2 



g(x)x n 3 dx . 
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which is an n x n linear system of equations for the coefficients Ck ■ In 
the code below the integrals on the right hand side are evaluated in 
terms of the function numericaLintegral. 

from numpy import arange, polyval, zeros, linalg 
f(x) = e~x 
interval = [-1, 1] 
nmax = 3 

data = zeros ((nmax, nmax)) 

rl , r2 = interval [0] , interval [1] 

for n in range(2, nmax+1) : 

a, b, c = zeros((n, n)), zeros(n), zeros(n) 
for j in ranged, n+1) : 

for k in ranged, n+1): 

a[j-l, k-1] = (r2~(2*n-k-j+l) - rl~ (2*n-k-j+l))/(2*n-k-j+l) 
b[j-l] = numerical_integral (f *x~ (n-j ) , rl, r2) [0] 
c = linalg. solve(a,b) 
h = (r2-rl)/40 
xs = arange (rl, r2+h, h) 
yl = [f(xi) for xi in xs] 
y2 = polyval (c, xs) 
err = abs(y2-yl) 
maxer = max (err) 

# Use trapezoidal rule to compute error 
intl = h*(sum(err) - 0.5*(err[0] + err[-l])) 

int2 = h*(sum(err~2) - . 5* (err [0] ~2 + err[-l]~2)) 

# Plots 

eplot = plot(f(x), (x, rl, r2) , color="black") 

polyplot = plot (polyval (c , x) , (x, rl, r2) , color="red", figsize=3) 
epoly = eplot + polyplot 

errplot = plot (abs (f (x) -polyval (c , x)), (x, rl, r2) , figsize=3) 

# Output text and graphics 
html("<hr>$n=7.s : $"7„n) 

html.table( [ ["$°/„s$"%latex(matrix(a) .n(digits=4) ) , 

"$°/,s$"7,latex(vector(b) .columnO .n(digits=4) ) , 
"$7,s$"7,latex(vector(c) .columnO .n(digits=4))]] , 
header=["$a$", "$b$", "$c$"]) 

html("$\\text{Abs . error = } °/ s\qquad\qquad\\text{L2 error = } 
%s$"% (maxer, int2)) 

html. table ( [["$\\text{Approximation (n = °/ s)}$"7 n, 

"$\\text{Abs. error function (n = %s)}$"7.n] , 



[epoly, errplot]], header=True) 



n = 2 : 



/ 0.5657 0.0000\ ( 0-7357A ( 1.104\ 
1^0.0000 2.00oj \ 2.35oj \l-17sj 



Abs. error = 0.439442311301 L2 error = 0.0530392393301 
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Approximation (n - 



Error functioning 2) 




0.4 












\ / 01 





re = 3 : 



^0.4000 0.0000 0.6667\ 
0-0000 0-8667 0-0000 ) 
^0.6667 0.0000 2.000/ 



O.S7S9\ /0.5367\ 
0.7357 1-104 ] 

2.350/ \0.9963/ 



Abs. errors 0.0816279606535 



L2 error = 0.0014675705349 



Approximation (n = 3) 



Error function(n = 3) 
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Figure 8. The picture shows the (n — l)th order poly- 
nomial approximation for the function e x on the interval 
[— 1, 1] in the cases of n = 2 and n = 3. 
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4. Concluding remarks 

During its initial years of development, the Sage project has grown 
to an environment which offers an attractive alternative for the com- 
mercial packages in several areas of computational mathematics. For 
the purpose of scientific computation teaching, the functionalities of 
Sage are practically the same as those of commercial packages. While 
free availability to instructional purposes is a very significant advan- 
tage, there are also other important factors from the learner's point of 
view: 

(1) The Python language can be used also for many other purposes 
not tied with the scientific computing. A wide selection of ex- 
tensions and other special libraries are available in the Internet. 

(2) The support of advanced data structures and support of object- 
oriented data types and modular program structure is available. 

(3) There is an active users' forum. 

It is likely that the Sage environment in education will become more 
popular on all levels of mathematics education from junior high school 
to graduate level teaching at universities. The support of symbolic 
computation via Maxima and various numerical packages are notewor- 
thy in this respect. For purposes of teaching scientific computing, the 
Sage environment and the modules it contains form an excellent option. 
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